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We construct efficient banks of templates suitable for all-sky narrow-band searches of almost 
monochromatic gravitational waves originating from spinning neutron stars in our Galaxy in data 
collected by interferometric detectors. We thus assume that both the position of the gravitational- 
wave source in the sky and the wave's frequency together with spindown parameters are unknown. 
In the construction we employ simplified model of the signal with constant amplitude and phase 
which is a linear function of unknown parameters. All our template banks enable usage of the fast 
Fourier transform algorithm in the computation of the maximum-likelihood J^-statistic for nodes 
of the grids defining the bank and fulfill an additional constraint needed to resample the data to 
barycentric time efficiently. Our template banks are suitable for larger range of search parameters 
than the banks previously proposed and compared to them they have smaller thicknesses for certain 
values of search parameters. 
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I. INTRODUCTION 

One of the primary sources of gravitational waves 
which currently operating ground-based interferomet- 
ric detectors LIGO JL], Virgo |2,], GEO600 Q, and 
TAMA300 Q are trying to detect, are rotating neutron 
stars located in our Galaxy. They are expected sources 
of almost monochromatic gravitational waves and in the 
present paper we consider the problem of construction of 
efhcient banks of templates needed to detect this kind of 
waves. 

Depending on what is a priori known about 
gravitational-wave sources we are looking for, searches 
can be splitted into targeted, directed, and all-sky (or 
blind) searches. In targeted searches both the position 
of the source in the sky and the wave's frequency to- 
gether with spindown parameters are known. If one as- 
sumes that only the position of the source in the sky 
is known but one does not know the frequency and the 
spindown parameters, one performs directed searches for 
gravitational-wave signals. Finally, in all-sky or blind 
searches one assumes that both the position of the source 
in the sky and the frequency and spindown parameters 
are not known. 

Several searches for almost monochromatic gravita- 
tional waves originating from spinning neutron stars in 
our Galaxy were already performed in the data collected 
by the LIGO, Virgo, and GEO600 detectors. The re- 
sults of targeted searches were published in Refs. 
(with recent searches for gravitational waves from Crab 
and Vela pulsars reported in and respectively). 
The results of the first directed search for gravitational 
waves from the supernova remnant Cassiopeia A were 
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published in Ref . TT] . Results of all-sky searches of data 
collected du ring LIGO science runs S2-S5 were reported 
in Refs. [la - flq . An Einstein@Home initiative (running 
on the BOINC infrastructure [13]) has performed all-sky 
searches on S4 and S5 LIGO data |18l - [20| . Recent results 
for the searches of continuous waves with the LIGO and 
Virgo detectors were shortly reviewed in (2ll |. 

In all these searches several different data analysis 
strategies were employed. In the present paper we re- 
strict ourselves to detection of gravitational-wave sig- 
nals and estimation of their parameters by means of 
the maximum-likelihood (ML) principle. We also assume 
that the noise in the detector is Gaussian and stationary 
(details of the ML detection in Gaussian noise can be 
found e.g. in [22,l23l] and in Chapter 6 of [2^]). In the case 
of targeted searches several statistics by means of which 
one can test whether data contains gravitational-wave 
signal were derived in Refs. [l^Hgj from the ML principle 
and also using the Bayesian approach together with the 
composite hypothesis testing. The ML princmle was also 
used in the directed search reported in Ref. (see also 
[131 )■ Data analysis tools and algorithms needed to per- 
form, within the ML approach, an all-sky narrow-band 
search for almost monochromatic gravitational- wave sig- 
nals were developed in detail in the series of papers |28l - 
m (see also Refs. [H,!!!). 

Within the ML approach one considers the likelihood 
ratio A[x; 6] which is a function of the data x and the pa- 
rameters 6 of the gravitational-wave signal we are look- 
ing for. ML detection of the signal relies on the com- 
putation of A[x; 6] maximized over all possible values of 
the parameters 6 and comparing this maximum with a 
threshold. In the case of directed or all-sky searches the 
unknown parameters 9 can be divided into two groups, 
6 = (A,^). The first group A consists of four extrinsic 
(or amplitude) parameters: an overall amplitude of the 
waveform, its initial phase, the polarization angle of the 
wave, and the inclination angle of the star's rotation axis 
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with respect to the hne of sight. The second group ^ con- 
tains intrinsic (or phase) parameters: the frequency of 
the wave, the spindown parameters, and the two more pa- 
rameters depending on the position of the gravitational- 
wave source in the sky (they are known in the case of 
directed searches). Maximization of the A with respect 
to amphtude parameters A can be done analyticahy by 
solving the set of equations dA/dA = 0; their solution 
with respect to A defines the ML estimators A — A[x; ^] 
of the parameters A. Then the J- -statistic is defined as 
the logarithm of the likelihood ratio A after replacing in 
A the amplitude parameters A by their ML estimators 
A: J^[x;i] :=lnA[a:;A[x;|],C]. 

The J^-statistic is a multidimensional random field as 
it depends on the detector's noise which is a stochastic 
process. Maximization of the J^-statistic over the phase 
parameters ^ can be done only numerically. To find the 
maximum of the J'-'-statistic one constructs a bank of tem- 
plates in the space of the parameters ^, which is deter- 
mined by a discrete set of points, i.e. a grid in the pa- 
rameter space. The grid is chosen in such a way, that 
for any possible gravitational- wave signal present in data 
there exists a grid point such that the expectation value 
of the J^-statistic computed for the parameters of this 
grid point is not less than a certain fixed minimal value. 

In the case of directed searches in the construction of 
the template banks one can use the simplified polyno- 
mial phase model of the gravitational-wave signal. This 
model was introduced in Ref. [l^ (see Sec. V D and Ap- 
pendix C there); in the model the signal's amplitude 
is constant and the phase is a polynomial function of 
time. It was found (in Sec. V E of [1^) that in the 
case of directed searches the polynomial phase model re- 
produces very well the covariance matrix (defined as the 
inverse of the Fisher matrix) for the ML estimators of the 
parameters of the exact gravitational-wave signal. The 
polynomial model was accepted in the search reported in 
Ref. fljj (see also [27j). In this search the phase of the 
gravitational-wave signal was modeled as a third-order- 
in-time polynomial (i.e. the frequency of the wave and 
its first and second spindown parameters were taken into 
account), and a template bank based on a body-centered 
cubic lattice was used. Efficient banks of templates for 
a second-order-in-time pol yno mial phase model was re- 
cently constructed in Ref. |35| . 

In the series of papers [28l - [3^ it was argued and 
checked by numerical simulations of the detection pro- 
cedure that in the case of all-sky narrow-band searches 
in the construction of the template banks one can employ 
another simplified model of the gravitational- wave signal, 
the so called linear phase model, in which the amplitude 
of the signal is assumed to be constant and the signal's 
phase is a linear function of the unknown parameters. 
The model was introduced in Sec. V B of Ref. [l^, and 
in Sec. V E of [1^ it was checked that the linear model 
reproduces well the covariance matrix (defined again as 
the inverse of the Fisher matrix) for the ML estimators 
of the parameters of the exact gravitational- wave signal. 



This model is accepted in the present paper. 

For the linear phase model the expectation value of the 
J^-statistic depends [see the key Eq. (|2.29p below] on the 
signal-to-noise ratio p and on the value of the autocovari- 
ance function Co(^,^') of the J^-statistic (the subscript 
'0' indicates that the autocovariance is calculated in the 
case when data is a pure noise) computed for the intrin- 
sic parameters of the template (^) and the gravitational- 
wave signal (^'), respectively. The signal-to-noise ratio 
p we can not control, therefore we fix its minimal value. 
Then to construct bank of templates one needs to choose 
some minimum value Cmin of the autocovariance func- 
tion Co and look for such a grid of points that for any 
point ^' in the intrinsic parameter space there exists a 
grid node ^ such that the autocovariance Co(^,^') com- 
puted for the parameters ^ and ^' is not less than Cmin- 
The autocovariance Cq{t) (for linear phase model it de- 
pends on ^, ^' only through the difference t :— ^ — 
can be approximated by taking the Taylor expansion (up 
to the second-order terms) of Cq around its maximum at 
T — 0. Then isoheights of such approximated autoco- 
variance function Cq are hyperellipsoids. The problem of 
constructing bank of templates can thus be formulated 
as a problem of finding optimal covering of the signal's 
parameter space by means of identical hyperellipsoids de- 
fined as isoheights of the autocovariance function Cq of 
the /"-statistic. 

In our paper we are interested in such searches for al- 
most monochromatic gravitational- wave signals for which 
the number of grid points in the parameter space is very 
large and the time needed to compute the J^-statistic for 
all grid nodes is long. Then it is crucial to use in the 
computation the fast numerical algorithms. Because the 
computation of the J^-statistic involves calculation of the 
Fourier transform, one would like to use the fast Fourier 
transform (FFT) algorithm. The FFT algorithm com- 
putes the values of the discrete Fourier transform (DFT) 
of a time series for a certain set of discrete frequencies 
called the Fourier frequencies. Thus it will be possible to 
use the FFT algorithm in computation of the J-"-statistic, 
provided the grid points will be arranged in such a way, 
that the frequency coordinates of these points will all co- 
incide with the Fourier frequencies. All grids constructed 
in our paper fulfill this requirement. 

The construction of efficient template banks for 
matched-filtering searches was considered in Ref. [s^ 
and the usage of random template banks and relaxed 
lattice coverings for gravitational-wave searches was re- 
cently discussed in Ref. 'st'I (see also [38-'40']). As ex- 
plained above we are interested in searches involving data 
streams so long, that the time performance of the search 
crucially depends on the ability of using the FFT al- 
gorithm. This enforces the above-mentioned constraint 
which is not fulfilled by grids considered in Refs. [36l - l40| . 
Therefore our work can be considered as beiiig comple- 
mentary to the studies performed in Refs. [36l - l40l |. Grids 
enabling the use of the FFT algorithm in the ML detec- 
tion of gravitational-wave signals from white-dwarf bi- 
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naries in the mock LISA data challenge were devised in 
Ref. dll (see also E^jlill), where the geometric approach 
(initialized in [il. l45j for searches of gravitational waves 
from inspiralling compact binaries and then developed 
also for searches of continuous gravitational waves, see 
[4^) was employed and the grids were constructed by 
some deformation of the optimal lattice coverings A*i in 
d = 3, 4 dimensions. The algorithm needed to construct 
bank of templates for all-sky narrow-band searches for 
almost monochromatic gravitational waves fulfilling the 
FFT-related constraint was devised in Sec. IV of Ref. 
[32j. In Sec. V below we compare the grids constructed 
in the present paper with those developed in {3^ . 

The organization of the paper is as follows. In Sec. II 
we introduce the linear phase model of the gravitational- 
wave signal. We consider here the phase with one spin- 
down parameter included. For this model we compute 
the J^-statistic and its expectation value in the case when 
the data contains the gravitational-wave signal. In Sec. 
Ill we introduce some mathematical notions related with 
coverings and we formulate the constraints we want to 
force on grids. Section IV is devoted to construction 
of two different families of grids. The grids enable us- 
age of the FFT algorithm in the computation of the J-- 
statistic for nodes of the grids and fulfill an additional 
constraint needed to resample the data to barycentric 
time efficiently. In Sec. V we discuss our results. Ap- 
pendix |X] gives some details of the optimal 4-dimensional 
lattice and Appendix [B] contains a brief sketch of the 
algorithm we use to find covering radius of given lattice. 

II. AUTOCOVARIANCE FUNCTION OF THE 
J^-STATISTIC 

We assume that the noise n in the detector is an ad- 
ditive, stationary, Gaussian, and zero-mean continuous 
stochastic process. Then the logarithm of the likelihood 
function is given by 

In A[x] = {x\h) - i {h\h), (2.1) 

where x denotes the data from the detector, h is the 
deterministic signal we are looking for in the data, and 
(• I • ) is the scalar product between waveforms defined by 

(/^,|M:=4Re^°°^i^g^d/. (2.2) 

Here ~ stands for the Fourier transform, * denotes com- 
plex conjugation, and Sn is the one-sided spectral density 
of the detector's noise n (Sn is defined thus for frequen- 
cies < / < +00). 

We are interested in narrow-band searches for almost 
monochromatic signals, i.e. such signals for which the 
modulus of the Fourier transform is concentrated (for 
frequencies / > 0) around some 'central' frequency fc > 
and Sn is a slowly changing function of / in the vicinity 



of the frequency fc- If both waveforms hi and /i2 in Eq- 
(|2.2p have their Fourier transforms concentrated around 
the same frequency fc, then we can replace Sn{f) in the 
integrand of (|2.2p by S'„(/c) and, after employing the 
Parseval's theorem, approximate the scalar product by 

ti+To/2 

(/iilM = ^;|^ / ^iW'^2(t)dt = ^^(/ii/z2). 

ti-To/2 

(2.3) 

Here (ij — To/2; ti + To/2) denotes observational interval, 
so To is the length of observation time, and ti — To/2 is 
the moment at which the observation begins. The time 
averaging operator ( ■ ) is defined by 

ti+To/2 

W:=^ J hit)dt. (2.4) 

ti-To/2 

Using the formula (|2.3p we can write the log likelihood 
ratio from Eq. (|2.ip as 

\nA[x]-^^^[{xh)-^{h^)). (2.5) 

In construction of template banks we employ an ap- 
proximate model of the continuous gravitational-wave 
signal from a rotating neutron star (this model was in- 
troduced in Sec. V B of Ref. [1^). The approximation 
relies on (i) assuming that the amplitude of the signal 
is constant, so we neglect the slowly varying modulation 
of the signal's amplitude due to motion of the detector 
with respect to the solar system barycenter (SSB); (ii) 
neglecting all spin downs in the phase modulation due 
to motion of the detector with respect to the SSB; (iii) 
discarding perpendicular to the ecliptic component of the 
vector connecting the SSB and the detector. This leads 
to the signal's model which is called linear because it has 
the property that its phase is a linear function of the pa- 
rameters. The approximate signal's model can be written 
as 

h{t; /lo, $0, €) = ho sin ($(<; |) + $0) , (2.6) 

where ho ia a constant amplitude and $0 is a constant 
initial phase. The time-dependent part ^{t; ^) of the 
phase depends on the s + 3 parameters ^, 

I = (wo,...,Ws,ai,a2), (2.7) 

and has the following form 

^) = I] ^fc (y) + "1 ^1 + "2 M2(i). (2.8) 

fe=0 ° 

The dimensioneless parameters ujk are defined as 
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where f^^^ = fo is an instantaneous frequency of the 

(k) 

gravitational wave computed at the SSB att — and /q 
(fc = 1, . . . , s) is the fcth time derivative of the instanta- 
neous gravitational- wave frequency at the SSB evaluated 
at t = 0. The parameters ai and a2 are related with 
the position of the gravitational-wave source in the sky 
through the definitions 

ai := 27r/o(sinQ!COS(5cos£ + sinJsine), (2.10a) 

q;2 := 27r/o cos a cos (5, (2.10b) 

where a is the right ascension and 6 is the declination of 
the source, e is the obliquity of the ecliptic. The functions 
fj-i(t) and fJ.2{t) are known functions of time, 

fii{t) :=i(i?|s(i) + Ry^{t) cose), (2.11a) 

/.2(i) :=i(i?lsW+^l'W), (2.11b) 

where (-RgS' ^eS' components of the vector 
joining the SSB with the center of the Earth in the SSB 
coordinate system, and {R^ , i?g , ) are the compo- 
nents of the vector joining the center of the Earth and the 
detector's location in the celestial coordinate system^ 
Let us introduce two new parameters. 



hi := /iocos$o, ^2 ft-osin^o- 



(2.12) 



Then Hq — ^Jh\ + h"^ and the signal h can be written as 
follows 

h{t;hi,h2,$) = /iisin$(t;|) + /i2COS$(t;|). (2.13) 
The time average {h"^) equals 

(h^) = ^{hl + hi) + ^{hj - /i?)(cos2$) + /ii/i2(sin2<I>). 

(2.14) 

For observations which last at least several hours and 
for gravitational- wave frequencies of the order of tens of 
Hertz or higher, to a good approximation 

(sin 2$) = 0, (cos 2$) = 0. (2.15) 

Then the time average (|2.14p simplifies to 



{h')-^ihl + hl). 



(2.16) 



Making use of (|2.16p we compute the optimal signal-to- 
noise ratio p for the signal (|2.13l) : 



piho)^^/im = J^^{h^)=ho.L^. (2.17) 



Snifc) 



Snifc)' 



^ The definitions of the SSB and celestial coordinate systems are 
given in Sec. II of Ref. [^. 



Substituting Eqs. and into (USD we get 

the following formula for the log likelihood ratio of the 
signal ((2l3t : 



lnA[a;;/ii,/i2,|] = 



2 To 

{fc. 



hi {x{t) sin$(t;|)) 



+ h2 {x{t) cos^t;^)) --(hl + hl) 

(2.18) 

Next we maximize In A with respect to the parameters 
hi and /i2 by solving equations 



d In A[a;; hi, /12, C] 
dhi 



0, i = l,2. 



The unique solution to these equations, 
/ii = 2(xsin$), /i2 = 2(xcos$) 



(2.19) 



(2.20) 



gives the maximum-likelihood estimators of the param- 
eters hi and /i2- After replacing in Eq. (|2.18p the pa- 
rameters hi and /12 by their estimators hi and /i2, we 
obtain the reduced log likelihood ratio which is called 
the T -statistic: 

I'lx;^] := lnA[a;; hi, /i2,l] 



Snifc) 



{x{t) sin'I>(t;^))^ + {x{t) cos$(i;|))^ 

(2.21) 



Making use of exp(— i"I>) — cos$ — isin$ (for $ e M) and 
the definition (|2.4p . it is easy to rewrite the J^-statistic 
in the form 



T[x;^] 



Sfi ifc)To 



ti+To/2 



J x(0exp(-i$(i;O 



ti-To/2 



dt 



(2.22) 



We now study the expectation value of the J^-statistic 
(|2.2ip in the case when the data x contains some 
gravitational- wave signal h, i.e. 



x{t) = n{t) + h{t;9'). 



(2.23) 



where 9' = {h'i,h'2,£,') collects the parameters of the 
gravitational-wave signal present in the data. We want 
thus to compute 

Ei{^[x; iW - E{ J-[n(t) + h{t; 6');^] }, (2.24) 

where the subscript '1' means that the expectation value 
is computed in the case when the data contains some 
signal. One can show that 

¥.i{F[x; i]} = l + \ p{h',)^ (( sin [$(i; |) - <J>(t; |')] )' 
+ (cos[<i>(t;0-<i>(i;O])'), (2.25) 
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where p(/io) is the signal-to-noise ratio from Eq. (|2.17p 
computed for the signal h{t]9') (so /iq = y/hJf+h^). 
The right-hand side of the above equation can be ex- 
pressed in terms of the autocovariance function Cq of the 
J^-statistic (the subscript '0' means here that the autoco- 
variance is computed in the case when the data contains 
only noise). In the signal- free case the J^-statistic T[n; $] 
is the random field which depends on the parameters ^, 
and its autocovariance function is defined as 



Coit E{[J-[n; ^] - mo(|)] [T[n; |'] - mo(|')]} , 

(2.26) 

where toq is the signal-free expectation value of 



mo(0 E{^[n;^]}. 



(2.27) 



In Sec. IV of Rcf. 3u| it was shown that the autocovari- 
ance function Co of the J^-statistic for the narrow-band 
gravitational-wave signal of the form (j2.13p can be ap- 
proximated by 

-f (cos[$(i;|)-$(f;^')])', (2.28) 
therefore the expectation value (|2.25p can be written as 



(2.29) 



The phase $ of the gravitational- wave signal (12.131) de- 
pends linearly on the parameters ^ [see Eq. (j2.8p ]. there- 
fore the autocovariance (I2.28P depends only on the dif- 
ferences between ^ and 

Co(|,r) = (sin$(t;|-0)' + (cos$(i;|-r))'. (2.30) 

If one introduces t := ^ — one can thus write 

Coir) ^ (cos$(i;T))2 + (sin $(t; t))2. (2.31) 

Let us note that Co attains its maximal value equal to 1 
for T = (i.e. for | = |'). 

We will further approximate the formula (12.311) for the 
autocovariance function in the case when \t\ <C 1. We 
will also restrict ourselves to the phase $ of the signal 
depending on one spindown parameter, so the phase is of 
the form [see Eq. ([27 



$(i;r) =wo^+cJif ^ ) +aifii{t) + a2fi2{t), (2.32) 



where r = (wo, wi, ai, a2). The approximation relies on 
expanding the right-hand side of Eq. (I2.3ip in Taylor se- 
ries around r = up to terms quadratic in t. Such com- 
puted autocovariance function we denote by Ca- Making 
use of the obvious equalities 



$(i;r = O)=0, 



dTkdn 



(t;x)=0, kj = l,...,4, 



we get 



4 4 



Co(r) ^ Ca(r) 1 - 5m ^fci 

k=l 1 = 1 



(2.33) 



(2.34) 



where T is the 4-dimensional reduced Fisher matrix with 
elements equal to 



' \dTkdTi/ \dTk'\dTi/' 



k,l = 1,...,4. 



(2.35) 

Let us note that because the phase $(i; x) is a linear 
function of the parameters r, the elements Tki of the 
Fisher matrix F are constant: they do not depend on the 
values of the parameters r. 

It is convenient to introduce the dimensionless quantity 
Xi and to replace the time t by the dimensionless variable 

X, 



t 



Xi rp 1 rp Xi 



(2.36) 



The observational interval (ij— To/2; ti+Tc,/2) of the time 
t is transformed, according to (|2.36p . into the interval 
(—1/2; 1/2) of unit length. Averaging with respect to 
the variable x is thus defined as 

1/2 

{{g{x))) J g{x)dx. (2.37) 

-1/2 

It is easy to see that for any function of time f{t) we have 

ti+To/2 

(/W> = ^ / f{t)dt 

ti-To/2 
1/2 

f{t{x))dx = {{f{t{xm, (2.38) 

1/2 



where [see Eq. t{x) = (x + Xi)To- Making use 

of the above introduced definitions, the Fisher matrix F 
with elements given in Eq. (|2.35p and for the phase $ 
defined in Eq. (|2.32p can be written as 



6 



r(xi) 



12 



180 



3 



((a^M2)) 



\ 



^32 
^42 



((^?)) - ((Mi))2 ((^lM2))-((m))((/^2)) 



(2.39) 



where 

f 23 = f32 - ((xVi)) + 2xi((x/ii)) - ^((mi)), (2.40a) 

f24 = f 42 - ((xV2» + 2Xi{{xfl2)) - ^((M2)). (2.40b) 

As we have indicated above, the elements of the reduced 
Fisher matrix T depend on the dimensionless parameter 
Xi (and on the time-dependent functions /ii and /i2 intro- 
duced in Eqs. (12. lip — they are determined by the motion 
of the detector with respect to the SSB). 



III. BANKS OF THE TEMPLATES 

To search for the gravitational-wave signal in detec- 
tor's noise we need to construct a bank of templates in 
the space of the parameters ^ on which the J^-statistic 
[given in Eq. (|2.2ip ] depends. The bank of templates is 
defined by a discrete set of points, i.e. a grid in the pa- 
rameter space chosen in such a way, that for any possible 
signal with parameters 6' = {h[, /12, ^') there exists a grid 
point ^ such that the expectation value of the J^-statistic, 
Ei{T[x;^]} ^E{T[n{t) + h{t;B');^]}, computed for the 
signal with parameters 9' and for the grid point ^, is not 
less than a certain fixed minimal value, assuming that the 
minimal vaule of the signal-to-noise p ratio is a priori 
fixed. From Eq. (j2.29p we see that this expectation value 
depends on the signal-to-noise ratio p and on the value 
C'o (€,€') of the noise autocovariance function computed 
for the intrinsic parameters ^ and ^' of the template and 
the gravitational-wave signal, respectively. 

To construct the bank of templates one thus needs to 
choose some minimum value Cmin of the autocovariance 
function Co and look for such a grid of points that for any 
signal with parameters ^' there exists a grid node ^ such 
that the autocovariance Co computed for the parameters 
^ and ^' is not less than Cmin, 

Co{t^')>Cmin. (3.1) 

We employ the linear model of the gravitational-wave 
signal for which the autocovariance Co(^,^') depends on 
^, ^' only through the difference ^ — , therefore (13. ip 
can be rewritten as 

Co{^-^')>C^in. (3.2) 



In the rest of this paper we will approximate the autoco- 
variance function Co by means of the formula (12.341) . i.e. 
we will use the approximate equality 

Co{i,i') = C,{ti')- (3.3) 

By virtue of Eq. (j2.34p the inequality p.2p can be written 
as 

4 

E f fc; iik - Q (6 - Ci) < 1 - Cmin, (3.4) 

k,l=l 

which for the fixed ^ is fulfilled by all points ^' which 
belong to an hyperellipsoid with the center located at ^ 
and with size determined by the value of Cmin- 

Wa want to find the optimal grid fulfilling the re- 
quirement (13.41) . i.e. the grid which consists of possibly 
smallest number of points. Thus the problem of find- 
ing the optimal grid is a kind of covering problem, i.e. 
the problem to cover the d-dimensional Euclidean space 
(or, in data analysis case, the bounded region of the 
space) by the smallest number of identical hypercllip- 
soids. The thorough exposition of the problem of cov- 
ering d-dimensional Euclidean space by identical hyper- 
spheres is given in Chap. 2 of Ref. j47j . 

We restrict ourselves to grids which are lattices, i.e. 
to grids with nodes which are linear combinations with 
integer coefficients of some basis vectors. If the vectors 
(Pi, . . . , Pd) are the basis vectors of a d-dimensional lat- 
tice, then a fundamental parallelotope is the subset of the 
R*^ consisting of the points 

AiPi + ... + AdPd, < Ai,...,Ad < 1. (3.5) 

A fundamental parallelotope is an example of a funda- 
mental region for the lattice, which when repeated many 
times fills the space with one lattice point in each copy. 

The quality of a covering can be expressed by the cov- 
ering thickness p which is defined as the average number 
of hyperellipsoids that contain a point in the space. For 
lattice coverings their thickness can be computed as 

volume of one hyperellipsoid 
volume of fundamental region 

Thickness of lattice covering of d-dimensional space M'' 
with identical hyperspheres of radius R reads 

^ dr(d/2)|detE|' ^ ' 
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where E is the matrix made of the basis vectors of the 
lattice. 

Another important notion is the covering radius of a 
lattice. Consider any discrete collection of points Q = 
{Qi, Q2, ■ • •} C R'*. The covering radius i? of Q is the 
least upper bound for the distance from any point x of 
M'' to the closest point Qi of the collection Q, 

R sup inf |x- Q,|. (3.8) 

Then identical hyperspheres of radius R centered at the 
points of Q will cover W^, and no hyperspheres of radius 
smaller than R will cover it. Around each point one 
defines its Voronoi cell, V{Qi), which consists of those 
point of R'' that are at least as close to as to any 
other Qj, 

y(Q,) := {xe R'^ : |x-Qi| < |x-Qj|for all j}. (3.9) 

The interiors of the Voronoi cells are disjoint. Each face 
of the Voronoi cell lies in the hyperplane midway between 
two neighboring points Q^. Voronoi cells are convex poly- 
topes whose union is the whole R'^. If the collection Q 
forms a lattice, then all the Voronoi cells are congruent. 

A. Constraints 

As we have already mentioned in Sec. I we are inter- 
ested in searches for gravitational- wave signals with very 
large number of grid points in the parameter space. Then 
the time needed to compute the J^-statistic for all grid 
nodes is long and we want to speed up this computa- 
tion by employing the FFT algorithm. As the FFT al- 
gorithm computes the values of the DFT of a time series 
for a certain set of discrete frequencies (the Fourier fre- 
quencies), it will be possible to use the FFT algorithm if 
the frequency coordinates of grid points will all coincide 
with the Fourier frequencies. We want to construct grids 
which fulfill this requirement. 

Le the data collected by a detector form a sequence of 
N samples 

Xu, u=l,...,N, (3.10) 

and let the sampling-in-time period be At. Then the 
DFT algorithm calculates the Fourier transform of the 
data with the frequency resolution A/ = l/{NAt). The 
resolution of the dimensionless frequency parameter ujq 
[introduced in Eq. (|2.9p ] is thus 

Aujo = 2nToAf ^ 2tt, (3.11) 

because TV At = To. 

It is possible to modify the DFT algorithm in such a 
way, that the frequency resolution p. lip changes. There 
exisit two types of such modifications (see Appendix in 
Ref. [m for more details) : (i) zero-padding of the data, 
which makes the DFT more dense; (ii) folding of the data. 



which diminishes the frequency resolution. If we add 
A'fft — A^ zeros to the A'' samples, so we make the Fourier 
transform of A^fft data points, then the frequency reso- 
lution is 

A^ 

Awo-27r- . (3.12) 

iVppx 

The case A^fft — 2N corresponds to usual zero-padding 
of the data, where for A^ data points we add N zeros. If 
we in turn fold A^-point data stream p times (so we finally 
have N/2P data points), then the frequency resolution 
changes to 

Awo==2Px27r, p=l,2,.... (3.13) 

Let (Pi, P2, P3, P4) be the basis vectors of a 4- 
dimensional lattice we consider. To use the FFT algo- 
rithm we need such grid that all nodes can be arranged 
along straight lines parallel to the ujq axis and the dis- 
tance between neighboring nodes along these lines must 
be equal to the frequency resolution of the FFT algo- 
rithm. We thus require that (say) the first basis vector 
of the grid we are looking for is of the forrr0 

Pi = (Ac^o,0,0,0)T, (3.14) 

where Aujq is the frequency resolution of the FFT algo- 
rithm we want to use [it is given in Eqs. p.l2p or (I3.13p ]. 

There is another important constraint to be met in all- 
sky searches. To reduce the computational time needed 
to resample the data to the so called barycentric time 
(see e.g. Sec. Ill D in [H, Sec. V A in [32], and [48^) we 
require that the 3rd and the 4th component of the (say) 
second basis vector vanish; 

P2 = (P21,P22,0,0)T. (3.15) 

B. Replacing hyperellipsoid-coverings by 
hypersphere-coverings 

It is convenient to replace the problem of finding the 
optimal covering of space by identical hyperellipsoids by 
the problem of finding the optimal covering of space by 
identical hyperspheres. Let us denote the original space 
of grid parameters by (fi C M^*, t G f2) and let the 
space of the transformed grid parameters be SI' (fi' C M.^, 
t' e O'). We are looking for the linear transformation 

t' = F(xi,au„)-T, (3.16) 

which transforms the hyperellipsoid of the constant value 
of the autocovariance function into the hypersphere of 



^ In the rest of the paper we treat all 4- vectors as column 4x1 
matrices. We will also use matrix notation with superscript "T" 
denoting matrix transposition and "■" denoting matrix multipli- 
cation. 
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unit radius, 
equal to Cn 
(1321)1 



This hyperellipsoid of the constant value 
n is determined by the equation [see Eq. 



• f(xi) • T - (1 - C,„i„) = 0, (3.17) 

whereas the equation of the unit hypersphere in the fl' 
space reads 

t'~^-t' -1 = 0. (3.18) 
After substituting ([3TT6)) into ((3l7|) we get 

. F(xi, C^inV ■ F(Xi, C^in) • T - 1 = 0. (3.19) 



Equations p.l9p and (I3.17P describe the same hyperel- 
lipsoid if and only if 

F(Xi, Cn^inV ■ F(Xi, C„,in) = p.^^ ,2 ^iXi), (3.20) 



where 



-R(Cmin) Vl — Cn 



(3.21) 



is the average radius of the hyperellipsoid (I3.17p . 

The Fisher matrix F is symmetric and [what can be 
shown by means of Eq. (|2.39p ] it is strictly positive def- 
inite, i.e. • F • T > for any t ^ 0. Therefore the 
equation (I3.20p can be interpreted as the Cholesky de- 
composition of the matrix T / E? , which states that there 
exists the unique upper triangular matrix F (with strictly 
positive diagonal elements) fulfilling Eq. (j3.20D . In the 
rest of the paper we will assume that the matrix F is the 
result of the Cholesky decomposition (so it is an upper 
triangular matrix). 

The elements of the upper triangular matrix F depend 
on the parameters Xi ^-iid Cmin and, through the evalu- 
ation of time averages needed to compute the elements 
of the Fisher matrix F [see Eq. (|2.39p ]. on the position 
vector of the detector with respect to the SSB during ob- 
servational interval. Making use of Eq. (I2.39P it is easy 
to show that the (1,1) element of the matrix F depends 
only on Cmin and it is equal to 



1 



2^3(1 - C„ 



(3.22) 



When the basis vectors of the grid in ft' space are 
found, one transform them into space by means of 
the inverse matrix F~^. If C is the generating matrix of 
the grid in fi' space, then the generating matrix C of the 
corresponding grid in fl space can be computed as 



C = C'.(F-i(xi,C„,i„))T. 



(3.23) 



Both matrices C and C are lower triangular (the matrix 
F""'^ is upper triangular so its transpose is lower triangular 
as well). 



IV. CONSTRUCTION OF GRIDS 

We will construct in the present section two families of 
grids which meet the two constraints (I3.14[) and (|3.15l) . 
Let us transform the first basis vector Pi of the grid we 
are looking for, which is fixed and given in Eq. p.l4p . 
from CI into f2' space, 



P'l := F(xi,C„,i„)-(Ac^o,0,0,0)T. 



(4.1) 



Because the matrix F is upper triangular, this leads to 

P'i = (Au;^,0,0,0)T (4.2) 

where, by means of Eq. (I3.22p . the length Alu'q of the 
basis vector P'^ in Q' space equals 



2^/3(1 - Cn 



(4.3) 



We also require that the second basis vector P2 of the 
grid in space fulfills the constraint p.lSp . Therefore in 
il' space the second basis vector has to have its 3rd and 
4th component equal to zero: 

P^ := FiXu Cn^in) ■ (^21,^22, 0, 0)T = (P2I, ^22, 0, 0)^. 

(4.4) 

All grids constructed in the present section in il' space 
fulfill the constraints and (lOl) . Constructions 

of grids depend [through Eq. (14.31) ] only on the values of 
Awq oraz Cmin- AH grids constructed in fl' space are 
built up from hyperspheres of unit radius. 

We will give below explicit numerical results for grids 
computed for some exemplary values of the parameters 
of the search. As the number of data points we take 



N = 344656. 



(4.5) 



This number corresponds to the observational interval of 
2 sidereal days sampled with the time period of 0.5 s. We 
consider the following lengths of the FFT: 



A^FFT — 2 



19 



(4.6) 



For these values of the search parameters the frequency 
resolution of the FFT according to Eq. (|3.12p reads, re- 
spectively. 



AujQ = 4.13044, 2.06522, 1.03261. 



(4.7) 



We will consider the minimum value of the autocovari- 
ance Cmin taken from the interval (0.70, 0.999), and 
sometimes as a reference value of Cmin we take Cmin — 
0.75. For this reference value of Cmin the length Awq of 
the basis vector P'^ computed by means of Eq. (|4.3p for 
the frequency resolution given in Eq. (14.61) equals, respec- 
tively. 



Auj'o = 2.38471, 1.19235, 0.596177. 



(4.8) 
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A. Optimal Al lattice 

Our constructions are based on some deformations of 
the 4-dimensional optimal lattice covering of space M.'^. 
This is so called A\ lattice [13], which is generated by 
the following matrix (made of the basis vectors Vi, . . . , V4 
with components arranged into rows of the matrix): 



V = 



/ 


V5 








\ 




2 


VIE 
2 










%/5 


VI 


/To" 







2 


2 


V 3 


V 


V5 






2V2 / 


2 


2~ 


2 



(4.9) 



This matrix generates covering of space by hypersheres 
of unit radii. Euclidean lengths of the basis vectors are: 



|vi| = |V2| = |V3| = V5, |V4|=V2. 

Thickness of the optimal lattice covering A\ is 
7r2 2^2 



2|detV| 5^5 



^ 1.765529. 



(4.10) 



(4.11) 



In the construction of the constrained grids given below 
the crucial role is played by the lengths of the lattice 
vectors for the optimal lattice A'^, i.e. the lengths of the 
vectors joining any two nodes of the lattice. The first 100 
smallest squares of the lengths, in ascending order, are 
listed in Appendix [X\ 



B. Grids Si for Cmin < C^in 

Our first construction will be valid (for the reasons 
explained at the end of the present subsection), for the 
given value of the frequency resolution Awq, only in the 
case when the minimum value Cmin of the autocovariance 
is less than C,„i„, where C*^^^^^ is defined in Eq. (I4.32p 
below. 

We begin the construction from replacing the basis vec- 
tors Va (a = 1, . . . , 4) of the lattice A\ by the new basis 
consisting of the following vectors: 



mi 

ma 
m4 



= Vl 
= V2 
= V3 
= V4. 



- V4, 

- V4, 

- V4, 



(4.12) 



Generating matrix of the lattice A4 with rows made of 
these vectors reads 



/m7\ 



mo 



Vmj; 



2 




2 



2\/2 ^ 
1 

2V2 
1 

2\/2 
1 

2\/2 



(4.13) 



Lengths of all vectors m^ are the same, 

|mi| = |m2| = Imal = |m4| = \/2. 



(4.14) 



Cosine of the angle and the angle itself between any two 
of the vectors ma equal 



cos;3 = -l/4, /3 = 1.82348. 



(4.15) 



Let us make the tails of the vectors ma to coincide with 
the origin of the coordinate system. Equations (I4.14[) and 
(I4.15|) imply then that the vectors ma coincide with the 
side edges of a 4-dimensional simplex whose base is a 
3-dimensional regular tetrahedron (with vertexes deter- 
mined by the heads of the vectors ma). The vector 



hi := ^(mi 



ma + m3 + m4) = 0,0,0,- 



1 



2V2 



(4.16) 

coincides with the height of this simplex. 

Grid points of the A^ lattice are determined by the 
vectors 

w(i, j, k, I) :~ i mi + j m.2 + k m.j, + I m4, (4-17) 

where ijj^kyl are any integers. We look for integers 
i,j,k,l such that the vector w{i,j,k,l) has length not 
less than Aui'q and as close to Awg as possible. There 
are many different choices leading to vectors of the same 
lengths. In this subsection we choose i=j = k = l = l 
what gives the following vector (parallel to the vector 
hi): 



ms = (0,0,0,-V2)^. 



(4.18) 



After replacing in the matrix (j4.13l) the vector m4 by 
ms, we get another matrix generating the lattice A^, 



2 



2^/2 ^ 
1 

2\/2 
1 

2\/2 

-V2 y 



(4.19) 



Now in three steps we will transform the matrix (j4.19p 
into a lower triangular matrix generating the same A*^ 
lattice. The height of the simplex built up from the basis 
vectors related to this new matrix will be parallel to the 
Wg axis. 

1. We exchange the order of the basis vectors: 1 O 4, 
2 -;-> 3. This leads to the matrix 



/ 





\ 

\ 2 







-V2 \ 
1 

2\/2 
1 

2\/2 
1 

2\/2 



(4.20) 
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2. We transform the basis vectors by means of the 
reflection with respect to the hyperplanes uj'q — cJ^ — 
[with normal vector (—1,0,0,1)^] andajj — a'^ — 
[with normal vector (0, —1, 1, 0)^]. This leads to the 
exchange of the columns: 1 O 4, 2 o 3, and gives 
the matrix 



/ -\/2 



V 



1 

2\/2 
1 

2\/2 
1 

2^2 



\ 







VI 

2 



(4.21) 



It is convenient to have generating matrix with all 
diagonal elements positive. To achieve this we re- 
flect the basis vectors of the matrix (j4.2ip with re- 
spect to the hyperplane Wg = 0. We get the gener- 
ating matrix N built up from the basis vectors ria 
(a = 1,2,3,4): 



n. 



"3 



/ \/2 



1 

2\/2 
1 

2\/2 
1 

2\/2 



\ 







2 



(4.22) 



Finally we deform the optimal lattice in the follow- 
ing way: we squeeze this lattice in the direction of the lo'q 
axis by the factor q < I, where 

q := Aw^/%/2. (4.23) 

The resulting grid is generated by the matrix 

Ci = N • Q, (4.24) 



where Q is the diagonal matrix with elements 



Q := 



(? 

10 

10 

1 



(4.25) 



The elements of the generating matrix Ci thus read 



Ci 



\ 



\ 4^^o — 2~ — 2~ T ' 



(4.26) 



Obviously the first two rows of this matrix form the vec- 
tors which meet the two constraints (j4.2p and (|4.4I) . Let 
us denote by Si the grid generated by the matrix Ci. 

The thickness of the grid 5*1 expressed in terms of Aoj'q 
is equal 



3/2 



By means of Eq. (|4.3p Auj'q can be expressed by search 
parameters Aujq and Cmin- Then the thickness of the 
grid Si can be written as 



6 2Vl^ 



As an example let us take Cn 



(4.28) 



0.75 and Aujq = 



2.06522. Then Auj'q = 1.19235, the coefficient q of Eq. 



equals q 

ments 



0.843122, and the matrix Ci has ele- 



Ci = 



/ 1.192353885 
0.2980884713 
0.2980884713 

V 0.2980884713 





1.369306394 
-0.4564354646 
-0.4564354646 





1.290994449 

-0.6454972244 1.118033989 



(4.29) 



whereas the thickness of the grid is 

PS, = 2.094038. (4.30) 

Equations (j4.23p and (14. 3p imply that the construction 
described above works only when Auj'q < \/2, i.e. when 
the resolution Aujq of the parameter ujq and Cmin fulfill 
the inequality (see Fig. [1]) 

AUJO < 2\/6v/l - Cmin. (4.31) 

Moreover, there is no need to deform the optimal lattice 
in the case when Auj'q — \/2, i.e. when the above 



inequality becomes the equality. For the fixed value of 
the parameter Aujq let us denote by C^^^^ the limiting 
value of the parameter Cmin- Then 

Aujo = 2V6V1 - C^^in or C;;^ = l-^A^^ (4.32) 

When Awo ^ 4.13044, 2.06522, 1.03261, this condition 
holds for C^-^ ^ 0.289146, 0.822287, 0.955572, respec- 
tively. If for the given value of the resolution Aujq one 
wants to consider Cmin > Cjjjjj^, then the contruction of 
the grid described in the next subsection can be used. 
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^min 

FIG. 1. The shaded region in the plane (Cmin, Acjo) is made of 
the points for which the construction of the grid Si described 
in Sec. IIV Bl is possible. Outside this region the construction 
from Sec. IIV Cl ean be used. 



C. Grids Si for C„,in > Ci^i^ 

In this subsection we will describe construction of grids 
Si valid also when Cmin > C^in- We start again from 
looking for such basis vector of the optimal lattice 
which has length as close to Awg as possible (but not 
less than Awq). Making use of the basis vectors of 
the lattice [they are rows of the matrix N from Eq. 
(I4.22P ]. grid points of the lattice can be written in the 
form 

q(i, j, k, /) := i ni + j n2 + fc ns + ? n4, (4.33) 

where i,j,kj are integers. We require that the chosen 
vector is such that one can take it and some three out of 
four basis vectors iIq to form a new basis of the lattice 
A^. Let us denote the chosen vector (with the length 
closest to Awq) by q and let its components be 

q= (a,6,c,d)"r. (4.34) 

We will rotate this vector to make it parallel to the Aojq 
axis. 

Let us first find any angle /3 satisfying equations 

b 



cos/3 



I , sin /3 = 



(4.35) 



and let us introduce the matrix Ri(/3) describing rotation 



in the 2-plane {uj[,a[) by the angle 
Ri(/3) : 



/ 1 

COS P — sin /3 

sin/3 cos/3 

Vo 1 



Then we define 

qi Ri (/3) • q = (a, 0, V^^ + c^, d f . 



(4.36) 



(4.37) 



Next we rotate the vector qi in the 2-plane {a[,a'2) by 
the angle 7 such that 



cos 7 



sin7 



VP 



To do this we introduce the rotation matrix 



d2 



R2(7) 



/ 1 

10 

COS 7 — sin 7 

\ sin 7 cos 7 



(4.38) 



(4.39) 



and define the vector 

q2 := R2(7) • qi - (a, 0, 0, Vb^ + + d2)T. (4.49) 

Finally we rotate the vector q2 in the 2-plane {ujQ,a2) 
employing the matrix 



R3(^) 



/cos (5 

1 

1 

V sin S 




(4.41) 



where the angle S is such that (here |q| 
Va^ + b'^ + c'^ + (P) 



cos (5 = 



sin 5 = — 



+ c2 + d2 



We define 



q3 := R3(<^)-q2 = (|q|,0,0,0)T. 



(4.42) 



(4.43) 



The vector qs is obviously parallel to the Wg axis. 

We now replace one of the basis vectors ria of the lattice 
A\ by the vector q. The vector q is a linear combination 
ia'n.a such that at least for one label b we have ib — 1. 
We replace nf, by the vector q and form a new basis for 
the lattice A\ made of q and the rest of the vectors ria. 
Then we apply the triple rotation Ri(/3) • R2(7) ■ ^^{5) 
to this new basis. Arranging components od the vectors 
into rows we get the generating matrix of the lattice A\ 
of the form 



No 



|q| 

Pi P2 P3 Pa 

n r2 ra r4 

Sl S2 S3 Si 



(4.44) 
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Making two further rotations in the 2-planes {uj[,a[) and 
{ai,a'2) we can make the ps and p4 components of the 
second basis vector to be zero. After this the generating 



matrix of the lattice A\ has the form 



/|q| 

P'l P2 

r[ r'j r'^ r'^ 

\ Si So Sa 



(4.45) 



Now we squeeze the optimal lattice in the direction 
of the lu'q axis. To do this we employ the diagonal matrix 



Qi = 



/ (71 

10 

10 

V 1 



where the squeezing factor 

qi AcJo/|q| 



(4.46) 



(4.47) 



is not greater than 1. The generating matrix of the grid 
reads 



Ci = Ni-Qi. 



(4.48) 



The first basis vector of this grid (which components form 
the first row of the matrix Ci) coincides with the con- 
straint (|4.2I) and the second one (with components taken 
from the second row of the matrix Ci) fulfills the con- 
straint (|44)) . 

For the chosen vector q the construction described 
above leads to covering only for Acjq < |q|. In the limit- 
ing case when Awq = |q| the grid 5*1 coincides with the 
optimal grid A\. The thickness of the grid 5*1 equals 



2^2 



2|detCi| 5^5 Aw^' 



(4.49) 



or, if one expresses Awq by search parameters Awq and 



473^2 |q|yr^ 



5\/5 



Awo 



(4.50) 



D. Grids S2 

The family of grids Si was obtained in the two pre- 
vious subsections as the result of squeezing the optimal 
4-dimensional lattice A^. Our starting point in the con- 
structions of grids Si was the lattice made of unit 
hypersheres, we have thus employed the lattice with 
covering radius equal exactly to 1. As consequence of 
squeezing the resulting grids 5*1 are made of unit hypersh- 
eres and have covering radius less than 1. In the present 
subsection we will introduce a one-parameter family of 
grids containing the grids 5*1 constructed for Cmin < C'min 



(or, equivalently, for Aujq < \/2) as a special case and ful- 
filling the constraints (j4.2p and (j4.4p . We expect that for 
some value of the parameter we will get grid with cov- 
ering radius equal exactly to 1, i.e. we will obtain grid 
with thickness smaller than the thickness of grid 5*1. We 
have found that the construction described below gives 
grids better than 5*1 not only when Auj'q < but also 
for some values of Auj'q greater than and close to v^. 

We start our construction from replacing the generat- 
ing matrix N of the lattice A^ [given in Eq. (|4.22p ] by the 
matrix built up from the following vectors: 



oi := ni n2 ns - n4, 
Oa ■■= Ha, a = 2,3,4. 



(4.51) 



The matrix with rows made of the components of the 
vectors Oa generates the lattice A^ and has the form 



= 



1 
1 

2V2 
1 

2V2 



2 






2 



(4.52) 



Let us make the tails of the vectors (|4.51l) to coincide with 
the origin of the coordinate system, then one can view 
on them as side edges of some simplex with the height 
determined by the vector 



if 



2V2 



0,0,0 



(4.53) 



Euclidean lengthes of the vectors Oa (a = 1,2,3,4) and 
the angles between any two of them are the same as for 
vectors nia (a = 1,2,3,4) introduced in (I4.12p [see Eqs. 
(|4.14p and (14.15^ ]. Let us introduce four vectors per- 
pendicular to the vector h2, 



aa := Oa h2, a = 1,2,3,4. 



(4.54) 



One easily checks that aa • h2 = (for a = 1, 2, 3, 4). We 
also introduce unit vectors 



|h2| 



|a„ 



a = 1,2,3,4. 



(4.55) 



We now deform the lattice A\ in the following way: In 
the simplex determined by the vectors Oq we replace the 
vectors Oq by the vectors bg (a = 1, 2, 3, 4) of the form 



ha{a, k) :— fccosah2-ffcsinaaa, 



1,2,3,4, (4.56) 



where A; > and < a < f (so that < cos a < 1) are 
some parameters. One easily checks that 



|bg(a,fc)| = fc, a = 1,2,3,4 
ba(a, k) ■ h2 



|ba(a,fc)||h2| 



cos a. 



1,2,3,4, 



(4.57a) 
(4.57b) 
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so k is the length of the vectors ha and a is the angle 
between any of the vectors ha and the vector h2 . We are 
thus changing the lengths of the side edges of the original 
A^-based simplex and we are also modifying the angle be- 
tween the side edges (and consequently the height of the 
simplex). Let us also define the following vector parallel 
to the Ldn axis: 



h.{a, k) :— k cos a h.2. 



(4.58) 



In the next step we replace one of the vectors b^, say 
bi, by the vector parallel to the uj'^ axis, this new vector 
we define as the sum of all four vectors ha- We thus get 
the following tetrad of vectors Ca : 



Ci(a, k) := Y.'^a=i ^a{a, k) = 4h(a, k), 
Cj(a, fc) := bQ(a, fc), a = 2, 3, 4. 



(4.59) 



The matrix built up from the vectors Cq (arranged into 
its rows) is lower diagonal. This matrix for k = \/2 and 
a = arccosl/4 ^ 1.31812 reproduces the matrix N from 
Eq. (j4.22p generating the optimal lattice A\. 

We require now that the length of the vector Ci fulfills 
the constraint 



|ci(a, fc)| = Aw^), 



(4.60) 



where Awq is given in Eq. (|4.3p . Because |ci(a, k)\ — 
4fccosa [see Eq. (|4.58p ]. from (|4.60p one can express the 
length fc as a function of the angle a: 



k{a) 



4 cos a 



(4.61) 



After substituting (j4.6ip into (|4.59p . the generating ma- 
trix of the grid built up from the vectors depends only 
on the angle a and can symbolically be written as 



C2(a) 



/ Ci{a,k{a)Y \ 
C2 (a, k{a)y 
C'i{a,k(a)Y 

\ Ci{a,k{a)y ) 



(4.62) 



The matrix C2(a) is lower diagonal. Let us note that the 
matrix C2(a) reproduces the matrix Ci (generating the 
grid 5*1 in the case Cmin < C'min) ^ora Eq. (|4.26p if one 
takes a = arctan('y30/Ati;g). The thickness of the grid 
generated by C2(a) equals 



P 



24v^7r^ cot^ a 



(4.63) 



Let us denote by 6*2(0?) the grid generated by the ma- 
trix C2(q!). We will now find the optimal value of the 



angle a, i.e. this value which minimizes the thickness of 
the grid 5*2 (a). The grid Si was obtained as the result 
of squeezing the optimal lattice A\ made of unit hyper- 
sheres, therefore its covering radius is less than 1. Let 



R 



1.5 



1.0 



0.5 



0.0 L 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 

a 

FIG. 2. Covering radius R of the grid 52 (a) as function of 
the angle a for Aa;o = 1.19235. 



us denote it by Rsi, Rsi < 1- We wiU find the value 
of the angle a for which the covering radius of the grid 
generated by the matrix €2(0:) will be equal to 1. 

The covering radius of the grid for the given value 
of the angle a we find by means of the algorithm 
sketched in Appendix [BJ This algorithm finds the vertex 
of the Voronoi cell inside the fundamental parallelotope 
spanned by the edge vectors (I4.56p . The covering radius 
is the distance of this vertex to the nearest vertex of the 
fundamental parallelotope. We have checked that the al- 
gorithm described in Appendix |B] works for angles a in 
the range (O, |) . 

The exemplary dependence of the covering radius of 
the grid 6*2(0;) on the angle a is depicted in Fig. [2] (for 
Awq = 1.19235). We see that the covering radius mono- 
tonically increases with the value of the angle a. One can 
numerically find the value Ofopt of the angle a for which 
the covering radius is equal to one, i?(aopt) = 1- In the 
case Aoj'n = 1.19235 this value reads 



aopt = 1.36530894. 



(4.64) 



The lengths of the simplex side edges computed [by 
means of Eq. (I4.6ip ] for this value of the angle a equal 



k{ao 



pt) 



1.46090061. The matrix C2(aopt) generating 



the best grid 6'2(aopt) reads 
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1.192353885 

, , I 0.298088471 1.430165676 0, , . 

L2(aoptj - I 0.298088471 -0.476721892 1.348373130 I ' ^ ' 

0.298088471 -0.476721892 -0.674186565 1.167725385, 



The thickness of the grid S'2(aopt) reads 
^2 



PS2(aopt) 



2|detC2(aopt)| 



1.83792. (4.66) 



It is better by ~12% than the thickness of the corre- 
sponding grid (computed for Auj'q = 1.19235), 



PSi - PS2(aopt) 
PSi 



^ 12.2%. 



(4.67) 



The thickness of the grid 5*2 is only ^4% larger than the 
thickness of the optimal lattice A^, 



PAX 



4.10%. 



(4.68) 



V. DISCUSSION 

We have constructed grids 5*1 and S2 and computed 
their thicknesses for different minimal values Cmin of 
the autocovariance function of the J^-statistic. We have 
taken the values of Cmin from the interval (0.70, 0.999) 
and we have made computations for the three different 
resolutions Awq of the dimensionless frequency parame- 
ter wq: Awo = 4.13044, 2.06522, 1.03261 [see Eqs. (|I3|) - 
(|4.8p and the text around them] . The thicknesses of the 
constructed grids are listed in Table U and are depicted in 
Fig- 13 Some slots in Table |T] which correspond to grids S2 
and high values of Cmin are empty — this is so because the 
construction of grids S2 described in Sec. lIVDl works, for 
the fixed value of Acjo, only up to some maximal value 

of Cmin- 

In Fig.|3]we have additionally indicated the thicknesses 
of the grids computed by means of the algorithm pre- 
sented in Sec. IV of Ref. [34I . Not for all combinations of 
Awo and Cmin we were able to construct grid using this 
algorithm; for some of them the algorithm gave no re- 
sult after a very long time (of the order of several hours) 
or just broke down. For the combinations (Aw^Cmin) 
for which we have constructed such grids, in almost all 
cases they have thickness equal to a good accuracy to 
the thickness of the corresponding Si grid constructed 
by us. This is by no means an obvious result, because 
the algorithm described in Sec. IV of Ref. [32,] is different 
from the algorithms devised in the present paper. We 
have in particular checked that the basis vectors of grids 
5*1 constructed by us are different from the correspond- 
ing basis vectors produced by the algorithm taken from 
Ref. [3"?] (obviously with the exception of the first basis 



vector which is fixed by the constraint imposed on the 
grids) . We can thus conclude that the family of grids Si 
devised by us is equivalent (in the sense of possessing the 
same thickness) to grids of Ref. [32j for these values of 
(AojOiCmin) for which the algorithm of Ref. [32I] works, 
and simultaneously it provides grids for those combina- 
tions of (AwojCmin) for which the algorithm of [32] is 
unable to produce grids. 

The reason of introducing, besides grids Si , the family 
of grids 5*2 is that for some values of search parameters 
grids ^2 have smaller thicknesses than grids 5*1. From 
Fig. 12] we see that the ranges of the parameter Cmin for 
which grids 6*2 have smaller thicknesses than grids 
are different for different frequency resolutions Awq. For 
the largest resolution we consider, Awq — 4.13044, the 
grids Si have thicknesses smaller than the grids iS'2 in 
the whole range 0.70 < Cmin < 0.999, so for this value 
of Awq there is no advantage in using grids S'2. In the 
worst case (for Cmin = 0.72) the thickness of the grid Si 
is 2.0730 (and it is ~17% larger than the thickness of the 
optimal A\ lattice). For the resolutions Acjq = 2.06522 
and 1.03261 the grids S'2 are in general better for smaller 
values of Cmin- In the case Awq = 2.06522 the worst case 
is for Cmin — 0.89, then the thickness of the grid S'2 is 
2.1226, ~20% more than the thickness of the A\ lattice. 
For the smallest considered resolution Awq = 1.03261 the 
thicker S2 grids correspond to smaller values of Cmin; the 
thickest one, for Cmin = 0.70, has the thickness 3.3074, 
~87% larger than the thickness of the A\ lattice. 

All grid constructions devised in Sec. IV of our pa- 
per depend on the search parameters Awq and Cmin 
only through the quantity Awq, see Eq. (j4.3l) . It means 
that any two grids (belonging to Si or to S2 family) de- 
picted on different panels of Fig. [3] which are determined 
by search parameters (Aa;g,Cj^i„) and (Awg, C^j^) such 

that At^i/yr^C^ = Ac^oV^/l - ^min- leads to the 
same value of Awg [see Eq. (|4.3I) ]. they thus both corre- 
spond to the same grid in U,' space. In Fig. 2] we have 
shown the dependence of the thicknesses of the grids Si 
and S'2 on the value of the parameter Awq . One can see 
that the thicknesses of the grids Si depicted in Fig. 2] split 
into several branches (better visible for smaller values of 
Awq). These branches correspond to different choices of 
the lattice vector q introduced in Sec. IV C. According 
to Eq. (j4.50p the thickness of the grid S'l is proportional 
to the length of q, so the sequence of branches visible in 
Fig. m is in fact the sequence of the lengths of possible 
lattice vectors for the optimal lattice A\ [the squares of 
these lengths are given in Eq. (lAip ]. The first branch 
(going from left to right) corresponds to [q] = -\/2 [see 
Eq. (jAll) ]. the second one to ]q] = -\/3, and so on. 
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FIG. 3. Covering thicknesses p of the grids Si (denoted by 
diamonds), 52 (circles), and the grids generated by means of 
the algorithm taken from Ref. [3^ (squares) as functions of 
Cmin for the three different frequency resolutions Awo of the 
search. The upper horizontal line corresponds to the thickness 
of the 4-dimensional hypercubic lattice (it equals = 4.9348) 
and the lower horizontal line denotes the thickness of the op- 
timal lattice AX (it equals = 1.7655). 



TABLE I. Covering thicknesses of the grids Si and S2 as 
functions of the minimal value Cmin of the autocovariance for 
the three different resolutions Aljq of the frequency parameter 
LJo- Aujo ^ 4.13044, 2.06522, 1.03261. 

Aluo ^ 4.13044 Acjo = 2.06522 Acjo = 1.03261 



C*min 






PS2 




PSi 




PS2 




PSi 




PS2 





70 


1 


8135 


4.5461 


2 


2939 


1 


9174 


4 


5878 


3 


3074 





71 


1 


7830 


4.4753 


2 


2553 


1 


9026 


4 


5107 


3 


2570 





72 


2 


0730 


4.4042 


2 


2161 


1 


8868 


4 


4322 


3 


2056 





73 


2 


0356 


4.3329 


2 


1762 


1 


8700 


4 


3524 


3 


1532 





74 


1 


9976 


4.2612 


2 


1355 


1 


8537 


4 


2710 


3 


1004 
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1 


9588 


4.1898 


2 


0940 


1 
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4 
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3 


0466 
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1 
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2 


0517 


1 


8229 


4 


1035 


2 


9918 





77 


1 


8788 


4.0459 


2 


0085 


1 


8089 


4 


0171 


2 


9361 





78 


1 


8375 


3.9740 


1 


9644 


1 


7961 


3 


9288 


2 


8792 





79 


1 


7953 


3.9019 


1 


9192 


1 
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3 
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2 
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80 


1 
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3 8303 


1 
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1 
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3 
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2 
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81 


1 
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1 
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1 
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3 
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2 


7009 





82 


1 
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1 


7769 


1 
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3 
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2 


6349 
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1 


9306 


3.6175 


2 
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1 
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3 
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2 
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1 
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2 
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1 
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3 


3505 


2 
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1 


8135 
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1 
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1 
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3 
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2 
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86 


1 
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1 
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1 


8230 


3 
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2 
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1 
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3 3535 


1 
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1 
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3 


0201 


2 


3041 





88 


1 
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3 2963 


1 
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1 


9638 


2 


9016 


2 


2450 





89 


1 


7707 


3.2457 


2 


1962 


2 
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2 


7781 


2 


1741 





90 


1 
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3.2054 


2 


0940 


2 


4537 


2 


6488 


2 


0989 





91 


1 


8315 


3.1811 


1 


9866 


3 


7346 


2 


5128 


2 


0295 





92 


1 


7769 


3.1342 


1 
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4 


6874 


2 


3691 


1 


9571 





93 


1 
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4 
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2 
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1 
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8135 
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3 
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2 
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1 
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1 
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6.7260 


1 
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3 
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1 
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1 
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1 
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1 


7769 


3 
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1 
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4 
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99 


1 
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1 
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3 
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1 
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3 
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1 
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1 
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4 
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1 


7769 


3 


4414 
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1 
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1 
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1 
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3 
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1 
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1 
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1 
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3 
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1 
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1 
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1 
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3 
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1 
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1 
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1 
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3 


1342 
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1 
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1 
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1 


7965 


2 


8233 
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1 
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1 


7916 






1 


7769 


3 


2085 





998 


1 


7670 




1 


7670 






1 


7769 


5 


5767 





999 


1 


7657 




1 


7694 






1 


7769 
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FIG. 4. Covering thickness p of the grids Si (diamonds) and S2 (circles) as functions of the quantity Al^q. We have shown 
here ah grids Si and S2 depicted in three panels of Fig. [3] as well as some more grids computed for smaller values of AajQ- The 
upper horizontal line corresponds to the thickness of the 4-dimensional hypercubic lattice (it equals = 4.9348) and the lower 
horizontal line denotes the thickness of the optimal lattice A\ (it equals = 1.7655). 



Appendix A: Lengths of the lattice vectors for the 
optimal lattice Al 

Let the origin of the coordinate system coincides with 
some node of the A\ lattice, then one can compute the 
lengths of the vectors joining the origin with other nodes 
of the lattice. The first 100 smallest squares of the 
lengths, in ascending order, equal 

2, 3, 5, 7, 8, 10, 12, 13, 15, 17, 15, 17, 18, 20, 22, 23, 25, 27, 28, 
30, 32, 33, 35, 37, 38, 40, 42, 43, 45, 47, 48, 50, 52, 53, 55, 57, 
58, 60, 62, 63, 65, 67, 68, 70, 72, 73, 75, 77, 78, 80, 82, 83, 85, 
87, 88, 90, 92, 93, 95, 97, 98, 100, 102, 103, 105, 107, 108, 
110, 112, 113, 115, 117, 118, 120, 122, 123, 125, 127, 128, 
130, 132, 133, 135, 137, 138, 140, 142, 143, 145, 147, 148, 
150,152,153,155,157,158,160,162,163,165,... . (Al) 



Appendix B: Algorithm finding covering radius 

The algorithm searches within fundamental parallelo- 
tope of given lattice vertexes of its Voronoi cell, i.e. the 



points which are the most distant from any vertex of the 
fundamental parallelotope. The algorithm makes use of 
the function, which for a given set of points lying in- 
side the fundamental parallelotope, computes the mini- 
mal distance of every point from all vertexes of the par- 
allelotope, and picks up this point for which the minimal 
distance achieves maximum. 

The algorithms works iteratively. Let the fundamental 
parallelotope of some c?-dimensional lattice be spanned 
by the vectors (vi, . . . , v^) and let (i = 1,2,.. .) be 
the position of the point picked up at the ith stage. The 
the points considered in the (z-l-l)th stage are determined 
by the formula 

d 

Hi = ni_i + ^ a„v™,i, a„ G {-2, 1,0, 1,2}, (Bl) 

m— 1 

where 

Vm.j:=^-^v„, i = l,2,.... (B2) 

The vector no initializing the algorithm coincides with 
the geometrical center of the fundamental parallelotope. 
Operation of the algorithm is illustrated in Fig. [5] for 
some 2-dimensional lattice. 
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vector no coincides with the center of the fundamental parallelotope. b) Around the head of the vector no we construct grid 
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stages of the algorithm. 
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